function y = power_mod_new_res(m, p_mn, coef_mn, coef_Im, CONSTS)

    a = CONSTS.a;
    k0 = CONSTS.k0;
    d = CONSTS.d;
    
    [~, ~, ~, ~, E_phi_m_border, Ez_in] = func_HE_inside(1, m, CONSTS, coef_mn, p_mn);
    N_mn = func_norm_res(m, CONSTS, coef_mn, p_mn);
    bessel_0 = besselj(0, k0*p_mn*d);

    y = -2*pi^4*a^2*real((abs(coef_Im))^2*bessel_0*conj(((E_phi_m_border^2)/(N_mn))*bessel_0));

end





